Finite element analysis is a powerful computer-aided tool for solving engineering and physical problems governed by partial differential equations. Finite element analysis is used in approximating any continuous physical characteristic of an object or a geometric region, such as temperature, heat, vibration, fluid flow, electric field, pressure, etc. Finite element analysis is specially important when the geometry of the object to be modeled is relatively complex, e.g. includes small details, since the differential equations for such applications become increasingly difficult to approximate accurately.
An initial step of a finite element analysis involves the construction of a mesh of finite elements covering the geometric region whose vertices are in turn points selected on the surface and in the interior of the object. Accurate and modern finite element analysis requires that a flexible mesh improvement/refinement/derefinement method, alternatively described as a point placement/elimination method, be used in order to produce a mesh that satisfies at least five criteria: (1) the overall point or vertex density is able to be specified and modified; (2) the density of vertices should increase and/or decrease in critical regions, (i.e., in those which have small features and concave corners or as needed in some time dependent problems); (3) badly-shaped finite elements should be avoided; (4) a mesh derived from a set of vertices should be able to be refined and/or derefined in certain areas by the user or the finite element analysis solving method; (5) the mesh should respect the object boundaries and interfaces.
It is known that the generation of a 3-D mesh suitable for finite element analysis has been one of the most time consuming steps in using the computer to analyze a complex engineering problem. Most previous work on 3-dimensional mesh improvement/refinement/derefinement and/or point placement/elimination has been heuristic in nature. For accurate analysis many approaches have required significant manual interaction between the user and the mesh generation system.
Several scientific articles considering mesh generation methods for use in finite element analysis have been published. In particular several papers appearing in the International Journal for Numerical Methods in Engineering consider the generation of a mesh by using the Delaunay triangulation. These papers are as follows: Cavendish et al., vol. 21, pp. 329-347, (1985); Frey, vol. 24, pp. 2183-2200, (1987); Schroeder et al., vol. 26, pp. 2503-2115, (1988); Schroeder et al. vol. 29, pp. 35-55, (1990); Golias et al., vol. 37, pp. 793-812, (1994); Field et al. vol. 31, pp. 413-425, (1991); Wheaterill et al., vol. 37, pp. 2005-2039 (1994);
Other papers considering Delaunay triangulation are the following: Baker, Engineering with Computers, vol. 5, pp. 161-175, (1989); Mitchell et al., IEEE Transactions on Magnetics, vol. 28, pp. 1751-1754 (1992); Ruppert, Journal of Algorithms, vol. 18, 548-585 (1995); Rebay, J. Comp. Physics, vol. 106, 125-138 (1993); Joe, Siam J. Sci. Comput. vol. 16, 1292-1307 (1995).
Each of the aforementioned papers employs a teaching of Delaunay to achieve planar triangulation and three-dimensional tetrahedrization. A tetrahedral mesh is defined as a "Delaunay tetrahedrization" if and only if, for each tetrahedron in the mesh, the sphere defined by the four vertices of the tetrahedron, called its circumsphere, contains no mesh vertices in its interior. Vertices on the circumsphere's boundary are permitted. The same property applies to two-dimensional surfaces: a surface triangulation is defined as a Delaunay triangulation if and only if, for each triangle in the mesh, a circle defined by the three vertices contains no other vertices in its interior.
However, it is well known to those skilled in the art that the ability of the 3-dimensional Delaunay technique to produce quality meshes strongly depends on the point placement method. Quality is a very important factor since a bad-quality mesh can result in inaccurate numerical approximation. Consequently bad-shaped tetrahedra are to be avoided, especially undesirable "slivers," that is tetrahedra formed by 4 almost coplanar vertices and whose longest-edge and smallest-edge are of "comparable size". Most of the aforementioned papers discuss some point placement strategies and some mesh improvement and/or mesh refinement tools.
Cavendish et al. (1985) implement such a triangulation by rejecting from the set of all possible triangles which might be formed, those with non-empty associated circles. Those triangles not rejected, form the Delaunay triangulation. Schroeder et al. (1988) apply those teachings to three-dimensional surfaces.
Frey (1987) teaches a method for selective refinement of an initial triangulation. Grading of the mesh is controlled by a node spacing function wherein a prospective node is inserted and its spacing from adjacent nodes is evaluated. Each new prospective node is also tested to see if its insertion would lead to a Delaunay triangulation with an acceptable degree of spacing at new node. Shroeder et al. (1990) use an octree method for producing a set of points representing the geometry of the object.
Golias et al. (1994) describe an element-based method to perform automatic mesh refinement producing almost-Delaunay triangulations as follows: the midpoint of the longest edge of each target element is chosen and all the elements having this edge in common are partitioned; the non-Delaunay mesh thus obtained is then subjected to the repeated application of a set of local topological transformations performed over all the elements of the mesh, followed by an overall node relaxation technique, until the equilibrium of the mesh is achieved.
Field et al. (1991) also use local transformations for producing almost-Delaunay meshes satisfying a local max-min solid angle criterion. Weatheril et al. (1994) select the points to be added to the mesh between the centroids of the elements of the mesh by using different point distribution functions interpolated to the interior vertices. Baker (1989) considers the addition of a new vertex near each sliver followed by a retriangulation step. Mitchell et al. (1992) use an element size function throughout the mesh and a local realignment of the nodes for improving the mesh. Joe (1995) also uses local topological transformations for producing improved almost Delaunay meshes.
Ruppert (1995) teaches a point insertion method which guarantees the construction of good-quality 2-D Delaunay triangulations by adding points in the circumcenter of the worst small-angled triangles of the current mesh. Rebay (1993) teaches an alternative point insertion method based on selecting a maximal non-accepted triangle which is also adjacent to an accepted triangle; then selecting a point situated over the segment that joins the circumcenter of the two triangles.
In the scientific literature considerable attention has been also paid to methods that consider the refinement of any general non-Delaunay mesh by performing selective longest-edge bisection of a set of appropriate elements of the mesh. In 2-dimensions, longest-edge bisection is obtained by partitioning the element by the line defined by the midpoint of the longest-edge and the opposite vertex. In 3-dimensions, longest-edge bisection is in turn obtained by partitioning the element by the plane defined by the midpoint of the longest-edge and the two opposite vertices. It has been shown that longest-edge refinement methods improve the point distribution by maintaining some small-angled triangles which depend on the quality of the initial mesh. The following papers consider the application of these ideas and the generation of a mesh therefor:
Rivara, International Journal for Numerical Methods in Engineering, vol. 20, pp. 745-756, (1984a); Rivara, ACM Trans. on Mathematical Software, vol. 10, pp. 242-264 (1984b); Rivara, In Accuracy Estimates and Adaptive Refinements in Finite Element Computations, John Wiley and Sons, Chichester, pp. 359-370 (1986); Rivara, International Journal for Numerical Methods in Engineering, vol. 24, pp. 1343-1354 (1987); Rivara, International Journal for Numerical Methods in Engineering, Vol. 28, pp. 2889-2906 (1989); Rivara et al, Communications on Applied Numerical Methods, vol. 8, pp. 281-290 (1992); Liu et al., SIAM J. Sci. Computing, vol. 16, pp. 1269-1291 (1995); Liu et al., Mathematics of Computation, vol. 65, pp. 1183-1200 (1996); Muthukrishnan,et al. AIAA Journal, vol. 33, pp. 928-932, (1995); Jones et al., Computing Systems in Engineering, vol. 5, pp. 297-309, 1994.
Rivara (1984a) teaches two kinds of longest-edge bisection refinement methods in 2-dimensions, which produce some neighbor's refinement and some intermediate non-valid meshes, and demonstrates their mathematical properties; Rivara (1984b and 1986) discusses the application of those teachings in adaptive multigrid finite element analysis; Rivara (1987) teaches the implementation of a 4-triangles method also based on the properties of longest-edge bisections which also produces some intermediate non-valid meshes; Rivara (1989) generalizes those teachings to the derefinement of nested triangulations based on assigning to each new vertex its level of creation and using this "label" recursively to identify a set of vertices, edges and triangles to be eliminated. The derefinement method produces some intermediate non-valid meshes.
Rivara et al. (1992) teach and discuss a three-dimensional longest-edge refinement method based on performing recursive bisection of the target elements and some of their largest neighbors. Liu et al. (1995 and 1996) consider two kinds of 3-dimensional methods based on bisection.
Rivara (1984b) and Rivara et al. (1992) teach mesh data structures based on associating to each vertex, respectively, the set of ordered neighbor vertices, and the set of the triangular neighbor faces.
Muthukrishnan et al. (1995) teach a variation of the longest-edge method of Rivara et al. (1992) based on performing a sort of the list of elements to be refined in ascending order of their longest-edges; associating a refinement ratio to these elements and distributing it to their descendants; inserting the adjacent elements of each selected element into the list at a position based on its longest edge; and using the order to perform the element division process.
Jones et al. (1994) teach a parallel scalable bisection method for the refinement of 2-dimensional triangulations based on the 4-triangles method of Rivara (1984a) under the assumption that the processors communicate through a common shared memory. Each target triangle to be refined is assigned to a processor, which keeps track of the neighbor's refinement in the mesh producing intermediate non-valid meshes. In order to avoid synchronization problems (since two different processors could create vertices at the same location when bisecting a triangle), the parallel method determines sequences of independent sets of triangles and refines the triangles in these sets in parallel.
The following recent papers consider the combination of longest-edge point insertion techniques and Delaunay triangulations: Rivara et al., International Journal for Numerical Methods in Engineering, vol. 40, 581-597 (1997a); Rivara, Proceedings 5th International Meshing Roundtable, Pittsburgh, Oct. 10-11 1996, pp. 77-86 (1996); Rivara, International Journal for Numerical Methods in Engineering, vol. 40, 3313-3324 (1997); Rivara et al., AMD-vol. 220, The American Society of Mechanical Engineers, pp. 1-8 (1997b).
Rivara et al. (1997a) consider a Delaunay triangulation in 2-dimensions and teach a method that for each target element perform the Delaunay insertion of a cluster of points defined by the longest-edge bisection method taught by Rivara (1984a). This method improves the point distribution and maintains some small-angled triangles which depend on the initial triangulation.
Rivara (1996 and 1997) teaches a new longest-edge point insertion method for 2-dimensional Delaunay meshes which always improves the triangulation. For each target triangle, the point to be inserted is found by following the sequence of neighbor increasing triangles until the last two greatest triangles in the sequence share a common longest-edge called a terminal-edge. Then the midpoint of this longest-edge is Delaunay inserted in the mesh. Rivara (1997) proves that the systematic use of this method over the small-angled triangles of the mesh, produces good-quality triangulations with smallest angles greater than or equal to 30.degree.. Rivara (1996 and 1997) also outlines a generalization of this method to 3-dimensions giving a rough 3-dimensional recursive method that, for each target element, follows a 3-dimensional propagation path associated with the longest-edge bisection method and performs the Delaunay insertion of a cluster of points over terminal-edges associated with the 3-dimensional propagation path. Rivara et al. (1997b) report numerical experimentation with both the 2-D and 3-D methods.
Rivara (1996 and 1997) also teaches recursive 2-dimensional and 3-dimensional longest-edge refinement methods for non-Delaunay quality meshes that produce the same meshes as the methods of Rivara (1984a) and Rivara et al. (1992) only partitioning sets of elements that share terminal edges of the mesh (all the elements of each set sharing the common terminal edge).
In U.S. Pat. No. 4,912,664 to Weiss et al., a 2-D approach to generating a mesh using Delaunay triangulation is described. Weiss et al. employ an expert-system based point placement method and an element quality criterion for point selection. In U.S. Pat. No. 4,933,889 to Meshkat et al., a further 2-D approach for generating a mesh using symmetric axis decomposition is disclosed. The method taught therein is not point based, but does employ the symmetric axis of a body to enable finite elements to be calculated.
In U.S. Pat. No. 5,214,752 to Meshkat et al., a 3-D automatic point placement method for generating a mesh using Delaunay triangulation is described. The method considers a subset of the vertices, and for each vertex selects the point to be inserted between points of intersection between a sphere of radius R centered at the vertex being considered, and the boundaries and edges. In addition, for eliminating sliver elements from the mesh, the center of the circumsphere corresponding to such element is inserted in the mesh.
In U.S. Pat. No. 5,315,537 to Blacker, an automatic quadrilateral surface discretization method and apparatus that produces mesh of all quadrilateral elements is disclosed. In U.S. Pat. No. 5,440,674 to Park, a 2-dimensional method for producing uniform and graded triangulations using 4 local operators is disclosed.
Further patents that disclose methods related with mesh generation are the following: Nackman, U.S. Pat. No. 4,797,842; Finnigan et. al., U.S. Pat. No. 5,345,490; Meshkat et al., U.S. Pat. 4,933,889; Shigyo et al., U.S. Pat. No. 4,941,114; Arakawa, U.S. Pat. No. 5,010,501; Arakawa, U.S. Pat. No. 5,398,307; Glassner, U.S. Pat. No. 5,428,717; Holmes, U.S. Pat. No. 5,497,451; Meshkat, U.S. Pat. No. 5,553,206; Akiyama, U.S. Pat. No. 5,774,696; Yamashita el al., U.S. Pat. No. 5,760,779; Yamamoto el al., U.S. Pat. No. 5,748,865; Strumulo el al., U.S. Pat. No. 5,729,670; Kumashiro, U.S. Pat. No. 5,677,846; Yokota, U.S. Pat. No. 5,617,332. Mesh generation can also be employed to enable the rendering of three-dimensional displays of various bodies, as described for example in: Falk, U.S. Pat. No. 4,888,713.
Despite the above, there is still much room for improvement to: improve mesh quality, refine/derefine a quality mesh, automatically control mesh grading and refinement/derefinement, automatically control point density by point insertion/elimination, eliminate undesirable slivers without using special techniques, and prevent badly-shaped mesh elements from appearing in the resulting mesh. It would be desirable to allow searching selected points to be inserted and/or eliminated considering only the distribution of increasing neighbor elements, not the Delaunay property of the mesh. It would desirable to allow the refinement/derefinement of the mesh avoiding the use of intermediate non-valid meshes. It would be desirable to allow the refinement/derefinement/improvement of the mesh by using non-recursive methods. It would be desirable to readily allow integrated improvement/refinement/derefinement of the mesh by taking full advantage of the properties of the geometric terminal-edge abstraction. It would be desirable to readily allow the parallel scalable refinement/derefinement of the mesh by locally modifying only sets of neighbor elements sharing a common edge (longest-edge) and requiring no interprocessor communication. It would be desirable to readily allow the management of Delaunay and non-Delaunay meshes. And it would be desirable to allow improved mesh data structures to implement integrated mesh improvement/refinement/derefinement techniques.